title: “infercnv_run” author: “Brian Haas” date: “9/5/2018” output: html_document: default pdf_document: default —
Create the InferCNV Object
infercnv_obj = CreateInfercnvObject(
raw_counts_matrix="sim.data",
annotations_file="sim.sample.annots.txt",
delim="\t",
gene_order_file="../example/full_dataset/gencode_v19_gene_pos.txt",
ref_group_names=c("normal") )Filtering genes
Removing those genes that are very lowly expressed or present in very few cells
# filter out low expressed genes
cutoff=1
infercnv_obj <- require_above_min_mean_expr_cutoff(infercnv_obj, cutoff)## INFO [2018-09-14 10:54:18] ::above_min_mean_expr_cutoff:Start
## INFO [2018-09-14 10:54:18] ::process_data:Averages (counts).
## INFO [2018-09-14 10:54:18] Removing 14262 genes from matrix as below mean expr threshold: 1# filter out bad cells
min_cells_per_gene=3
infercnv_obj <- require_above_min_cells_ref(infercnv_obj, min_cells_per_gene=min_cells_per_gene)## INFO [2018-09-14 10:54:18] no genes removed due to min cells/gene filter## for safe keeping
infercnv_orig_filtered = infercnv_obj
#plot_mean_chr_expr_lineplot(infercnv_obj)
save('infercnv_obj', file = 'infercnv_obj.orig_filtered')Normalize each cell’s counts for sequencing depth
infercnv_obj <- infercnv:::normalize_counts_by_seq_depth(infercnv_obj)perform Anscombe normalization
Suggested by Matan for removing noisy variation at low counts
infercnv_obj <- infercnv:::anscombe_transform(infercnv_obj)
save('infercnv_obj', file='infercnv_obj.anscombe')log transform the normalized counts:
infercnv_obj <- log2xplus1(infercnv_obj)## INFO [2018-09-14 10:54:19] transforming log2xplus1()save('infercnv_obj', file='infercnv_obj.log_transformed')Apply maximum bounds to the expression data to reduce outlier effects
threshold = mean(abs(get_average_bounds(infercnv_obj)))
infercnv_obj <- apply_max_threshold_bounds(infercnv_obj, threshold=threshold)## INFO [2018-09-14 10:54:20] ::process_data:setting max centered expr, threshold set to: +/-: 4.1252634394848Initial view, before inferCNV operations:
plot_cnv(infercnv_obj,
output_filename='infercnv.logtransf',
x.range="auto",
title = "Before InferCNV (filtered & log2 transformed)",
color_safe_pal = FALSE,
x.center = mean(infercnv_obj@expr.data))knitr::include_graphics("infercnv.logtransf.png")
perform smoothing across chromosomes
infercnv_obj = smooth_by_chromosome(infercnv_obj, window_length=101, smooth_ends=TRUE)## INFO [2018-09-14 10:54:47] ::smooth_window:Start.save('infercnv_obj', file='infercnv_obj.smooth_by_chr')
# re-center each cell
infercnv_obj <- center_cell_expr_across_chromosome(infercnv_obj, method = "median")## INFO [2018-09-14 10:54:49] ::center_smooth across chromosomes per cellsave('infercnv_obj', file='infercnv_obj.cells_recentered')plot_cnv(infercnv_obj,
output_filename='infercnv.chr_smoothed',
x.range="auto",
title = "chr smoothed and cells re-centered",
color_safe_pal = FALSE)knitr::include_graphics("infercnv.chr_smoothed.png")
subtract the reference values from observations, now have log(fold change) values
infercnv_obj <- subtract_ref_expr_from_obs(infercnv_obj, inv_log=TRUE)## INFO [2018-09-14 10:55:19] ::subtract_ref_expr_from_obs:Start
## INFO [2018-09-14 10:55:19] inverting log2xplus1()
## INFO [2018-09-14 10:55:19] subtracting mean(normal) per gene per cell across all data
## INFO [2018-09-14 10:55:59] transforming log2xplus1()save('infercnv_obj', file='infercnv_obj.ref_subtracted')plot_cnv(infercnv_obj,
output_filename='infercnv.ref_subtracted',
x.range="auto",
title="ref subtracted",
color_safe_pal = FALSE)knitr::include_graphics("infercnv.ref_subtracted.png")
invert log values
Converting the log(FC) values to regular fold change values, centered at 1 (no fold change)
This is important because we want (1/2)x to be symmetrical to 1.5x, representing loss/gain of one chromosome region.
infercnv_obj <- invert_log2(infercnv_obj)## INFO [2018-09-14 10:56:29] invert_log2(), computing 2^xsave('infercnv_obj', file='infercnv_obj.inverted_log')plot_cnv(infercnv_obj,
output_filename='infercnv.inverted',
color_safe_pal = FALSE,
x.range="auto",
x.center=1,
title = "inverted log FC to FC")knitr::include_graphics("infercnv.inverted.png")
Removing noise
infercnv_obj <- clear_noise_via_ref_mean_sd(infercnv_obj, sd_amplifier = 1)## INFO [2018-09-14 10:56:57] :: **** clear_noise_via_ref_quantiles **** : removing noise between bounds: 0.938873836066226 - 1.06112616393377save('infercnv_obj', file='infercnv_obj.denoised')plot_cnv(infercnv_obj,
output_filename='infercnv.denoised',
x.range="auto",
x.center=1,
title="denoised",
color_safe_pal = FALSE)knitr::include_graphics("infercnv.denoised.png")
Remove outlier data points
This generally improves on the visualization
infercnv_obj = remove_outliers_norm(infercnv_obj)## INFO [2018-09-14 10:57:24] ::remove_outlier_norm:Start out_method: average_bound lower_bound: NA upper_bound: NA
## INFO [2018-09-14 10:57:24] ::remove_outlier_norm using method: average_bound for defining outliers.save('infercnv_obj', file="infercnv_obj.outliers_removed")plot_cnv(infercnv_obj,
output_filename='infercnv.outliers_removed',
color_safe_pal = FALSE,
x.range="auto",
x.center=1,
title = "outliers removed")knitr::include_graphics("infercnv.outliers_removed.png")
Find DE genes by comparing the mutant types to normal types, BASIC
Runs a t-Test comparing tumor/normal for each patient and normal sample, and masks out those genes that are not significantly DE.
load('infercnv_obj.final')
plot_data = infercnv_obj@expr.data
high_threshold = max(abs(quantile(plot_data[plot_data != 0], c(0.05, 0.95))))
low_threshold = -1 * high_threshold
infercnv_obj2 <- infercnv:::mask_non_DE_genes_basic(infercnv_obj, test.use = 't', center_val=1)## INFO [2018-09-14 10:57:50] Finding DE genes between Apt25 and normal
## INFO [2018-09-14 10:57:51] Found 605 genes as DE
## INFO [2018-09-14 10:57:51] Finding DE genes between Bpt5 and normal
## INFO [2018-09-14 10:57:53] Found 554 genes as DE
## INFO [2018-09-14 10:57:53] Finding DE genes between C1pt5 and normal
## INFO [2018-09-14 10:57:54] Found 502 genes as DE
## INFO [2018-09-14 10:57:54] Finding DE genes between D1pt75 and normal
## INFO [2018-09-14 10:57:55] Found 574 genes as DEplot_cnv(infercnv_obj2,
output_filename='infercnv.non-DE-genes-masked',
color_safe_pal = FALSE,
x.range=c(low_threshold, high_threshold),
x.center=1,
title = "non-DE-genes-masked")knitr::include_graphics("infercnv.non-DE-genes-masked.png")